# Load packages
library(tidyverse) # for data manipulation
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## Warning: package 'tibble' was built under R version 4.0.2
## Warning: package 'dplyr' was built under R version 4.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(rstanarm) # For bayesian estimates of alpha and beta diversity
## Loading required package: Rcpp
## This is rstanarm version 2.21.2
## - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
## - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
## - For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores())
library(gridExtra) # For arranging ggplots
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
library(mgcv) # For beta distribution (beta diversity)
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.
library(RColorBrewer) # colors for barplots
library(grid) # for text grobs in gridExtra
library(vegan) # for adonis
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.5-6
library(MASS) # for fitdistr
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(car) # for Anova
## Loading required package: carData
## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4
##
## Attaching package: 'car'
## The following object is masked from 'package:rstanarm':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
library(betareg) # for beta regression
library(lme4) # for lmer
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
##
## lmList
source("msc_codebits_Baysian.R")
# Read in mapping files and OTU tables
load("./3_5sp_mapping_otu_downstream/mf_alt_filt_final.RData")
load("./3_5sp_mapping_otu_downstream/otu_filt.RData")
load("./3_5sp_mapping_otu_downstream/otu_filt_inhibOnly.RData")
load("./3_5sp_mapping_otu_downstream/braycurtis_filt.Rdata")
load("./3_5sp_mapping_otu_downstream/unweighted_unifrac_filt.Rdata")
load("./3_5sp_mapping_otu_downstream/weighted_unifrac_filt.Rdata")
load("./3_5sp_mapping_otu_downstream/alpha_metrics.Rdata")
load("./3_5sp_mapping_otu_downstream/beta_metrics.Rdata")
mf_con <- mf_alt_filt_final %>%
filter(Bd_exposure == "Control")
mf_treat <- mf_alt_filt_final %>%
filter(Bd_exposure == "Bd-exposed")
dir.create("./4_Bayesian_models")
## Warning in dir.create("./4_Bayesian_models"): './4_Bayesian_models' already
## exists
#### PLOTTING EXP DESIGN ####
options(repr.plot.width=5, repr.plot.height=10)
mf_alt_filt_final %>%
mutate(Contaminated = factor(ifelse(orig_contam ==1, "!Contaminated",NA), levels=c("!Contaminated"))
, Bd_logload = (Bd_load)) %>%
ggplot(aes(x=time, y=indiv)) +
geom_line(aes(group=indivID, col=Bd_exposure)) +
geom_point(aes(group=indivID,bg=Bd_logload), cex=4, pch=21)+
scale_color_manual(values=c("black","blue","orange")) +
scale_fill_gradient(low = "white", high = "red") +
geom_vline(aes(xintercept=5.5), col="orange")+
geom_point(aes(group=indivID, col=Contaminated), cex=1, pch=19)+ ## NEW LINE
facet_wrap(~species, nrow=5) +
xlab("Time") +
ylab("Individual Toad")
## Warning: Removed 515 rows containing missing values (geom_point).
#### PLOTTING BETA COMPOSITION PLOTS ####
# What do different species look like? (Control only)
options(repr.plot.width=12, repr.plot.height=12)
gg_control_bc_sp <- mf_con %>%
ggplot(aes(x=NMDS1_braycurtis, y=NMDS2_braycurtis)) +
geom_point(aes(col=species), cex=3) +
ggtitle(label = paste0("Bray-Curtis (Stress:",as.character(round(mean(mf_alt_filt_final$NMDS_stress_braycurtis),0)/100),")"))
gg_control_uwu_sp <- mf_con %>%
ggplot(aes(x=NMDS1_unweighted_unifrac, y=NMDS2_unweighted_unifrac)) +
geom_point(aes(col=species), cex=3) +
ggtitle(label = paste0("Unweighted Unifrac (Stress:",as.character(round(mean(mf_alt_filt_final$NMDS_stress_unweighted_unifrac),0)/100),")"))
gg_control_wu_sp <- mf_con %>%
ggplot(aes(x=NMDS1_weighted_unifrac, y=NMDS2_weighted_unifrac)) +
geom_point(aes(col=species), cex=3) +
ggtitle(label = paste0("Weighted Unifrac (Stress:",as.character(round(mean(mf_alt_filt_final$NMDS_stress_weighted_unifrac),0)/100),")"))
gg_control_bc_time <- mf_con %>%
ggplot(aes(x=NMDS1_braycurtis, y=NMDS2_braycurtis)) +
geom_point(aes(col=time), cex=3)
gg_control_uwu_time <- mf_con %>%
ggplot(aes(x=NMDS1_unweighted_unifrac, y=NMDS2_unweighted_unifrac)) +
geom_point(aes(col=time), cex=3)
gg_control_wu_time <- mf_con %>%
ggplot(aes(x=NMDS1_weighted_unifrac, y=NMDS2_weighted_unifrac)) +
geom_point(aes(col=time), cex=3)
grid.arrange( gg_control_bc_sp,gg_control_bc_time
, gg_control_uwu_sp,gg_control_uwu_time
, gg_control_wu_sp,gg_control_wu_time
, nrow=3)
options(repr.plot.height=12, repr.plot.width=6)
gg_all_bc <- mf_alt_filt_final %>%
ggplot(aes(x=NMDS1_braycurtis, y=NMDS2_braycurtis)) +
geom_point(aes(col=species,pch=Bd_exposure ), cex=2) +
ggtitle(label = paste0("Bray-Curtis (Stress:",as.character(round(mean(mf_alt_filt_final$NMDS_stress_braycurtis),0)/100),")"))
gg_all_uwu <- mf_alt_filt_final %>%
ggplot(aes(x=NMDS1_unweighted_unifrac, y=NMDS2_unweighted_unifrac)) +
geom_point(aes(col=species,pch=Bd_exposure ), cex=2) +
ggtitle(label = paste0("Unweighted Unifrac (Stress:",as.character(round(mean(mf_alt_filt_final$NMDS_stress_unweighted_unifrac),0)/100),")"))
gg_all_wu <- mf_alt_filt_final %>%
ggplot(aes(x=NMDS1_weighted_unifrac, y=NMDS2_weighted_unifrac)) +
geom_point(aes(col=species,pch=Bd_exposure ), cex=2) +
ggtitle(label = paste0("Weighted Unifrac (Stress:",as.character(round(mean(mf_alt_filt_final$NMDS_stress_weighted_unifrac),0)/100),")"))
grid.arrange(gg_all_bc,gg_all_uwu, gg_all_wu
, nrow=3)
# # Check if marginal has an effect
# adonis2(braycurtis_filt ~ species:time, data=mf_alt_filt_final, by="margin")
# # Check main effects
# adonis2(braycurtis_filt ~ species + time + species:time, data=mf_alt_filt_final)
# # Include Bd_exposure
# adonis2(braycurtis_filt ~ species*time*Bd_exposure, data=mf_alt_filt_final)
#
# # Check if marginal has an effect
# adonis2(unweighted_unifrac_filt ~ species:time, data=mf_alt_filt_final, by="margin")
# # Check main effects
# adonis2(unweighted_unifrac_filt ~ species + time + species:time, data=mf_alt_filt_final)
# # Include Bd_exposure
# adonis2(unweighted_unifrac_filt ~ species*time*Bd_exposure, data=mf_alt_filt_final)
#
# # Check if marginal has an effect
# adonis2(weighted_unifrac_filt ~ species:time, data=mf_alt_filt_final, by="margin")
# # Check main effects
# adonis2(weighted_unifrac_filt ~ species + time + species:time, data=mf_alt_filt_final)
# # Include Bd_exposure
# adonis2(weighted_unifrac_filt ~ species*time*Bd_exposure, data=mf_alt_filt_final)
options(repr.plot.height=5, repr.plot.width=12)
# Bray curtis
mf_alt_filt_final %>%
mutate(Exposure=factor(prepost, levels=c("Pre","Post"))) %>%
ggplot(aes(x=NMDS1_braycurtis,y=NMDS2_braycurtis)) +
geom_point(aes(bg=Bd_exposure, col=Exposure), cex=2, alpha=0.8, pch=21) +
facet_grid(Bd_exposure ~ species) +
scale_color_manual(values=c("white","black")) +
scale_fill_manual(values=c("blue","red"))
# Unweighted Unifrac
mf_alt_filt_final %>%
mutate(Exposure=factor(prepost, levels=c("Pre","Post"))) %>%
ggplot(aes(x=NMDS1_unweighted_unifrac,y=NMDS2_unweighted_unifrac)) +
geom_point(aes(bg=Bd_exposure, col=Exposure), cex=2, alpha=0.8, pch=21) +
facet_grid(Bd_exposure ~ species) +
scale_color_manual(values=c("white","black")) +
scale_fill_manual(values=c("blue","red"))
# Weighted Unifrac
mf_alt_filt_final %>%
mutate(Exposure=factor(prepost, levels=c("Pre","Post"))) %>%
ggplot(aes(x=NMDS1_weighted_unifrac,y=NMDS2_weighted_unifrac)) +
geom_point(aes(bg=Bd_exposure, col=Exposure), cex=2, alpha=0.8, pch=21) +
facet_grid(Bd_exposure ~ species) +
scale_color_manual(values=c("white","black")) +
scale_fill_manual(values=c("blue","red"))
options(repr.plot.height=5, repr.plot.width=12)
# Bray curtis
mf_alt_filt_final %>%
mutate(Infection_status=factor(ifelse(PABD==0,"Not infected","Infected"), levels=c("Not infected","Infected"))) %>%
ggplot(aes(x=NMDS1_braycurtis,y=NMDS2_braycurtis)) +
geom_point(aes(bg=Bd_exposure, col=Infection_status), cex=2, alpha=0.8, pch=21) +
facet_grid(Bd_exposure ~ species) +
scale_color_manual(values=c("white","black")) +
scale_fill_manual(values=c("blue","red"))
# Unweighted Unifrac
mf_alt_filt_final %>%
mutate(Infection_status=factor(ifelse(PABD==0,"Not infected","Infected"), levels=c("Not infected","Infected"))) %>%
ggplot(aes(x=NMDS1_unweighted_unifrac,y=NMDS2_unweighted_unifrac)) +
geom_point(aes(bg=Bd_exposure, col=Infection_status), cex=2, alpha=0.8, pch=21) +
facet_grid(Bd_exposure ~ species) +
scale_color_manual(values=c("white","black")) +
scale_fill_manual(values=c("blue","red"))
# Weighted Unifrac
mf_alt_filt_final %>%
mutate(Infection_status=factor(ifelse(PABD==0,"Not infected","Infected"), levels=c("Not infected","Infected"))) %>%
ggplot(aes(x=NMDS1_weighted_unifrac,y=NMDS2_weighted_unifrac)) +
geom_point(aes(bg=Bd_exposure, col=Infection_status), cex=2, alpha=0.8, pch=21) +
facet_grid(Bd_exposure ~ species) +
scale_color_manual(values=c("white","black")) +
scale_fill_manual(values=c("blue","red"))
options(repr.plot.height=3, repr.plot.width=12)
# Bray curtis
mf_alt_filt_final %>%
mutate(Infection_status=factor(ifelse(PABD==0,"Not infected","Infected"), levels=c("Not infected","Infected"))) %>%
filter(prepost=="Post") %>%
ggplot(aes(x=NMDS1_braycurtis, y=NMDS2_braycurtis)) +
geom_point(aes(bg=time, col=Infection_status), cex=3, pch=21) +
scale_color_manual(values=c("white","red"))+
facet_wrap(~species, nrow=1)
# Unweighted Unifrac
mf_alt_filt_final %>%
mutate(Infection_status=factor(ifelse(PABD==0,"Not infected","Infected"), levels=c("Not infected","Infected"))) %>%
filter(prepost=="Post") %>%
ggplot(aes(x=NMDS1_unweighted_unifrac, y=NMDS2_unweighted_unifrac)) +
geom_point(aes(bg=time, col=Infection_status), cex=3, pch=21) +
scale_color_manual(values=c("white","red"))+
facet_wrap(~species, nrow=1)
# Weighted Unifrac
mf_alt_filt_final %>%
mutate(Infection_status=factor(ifelse(PABD==0,"Not infected","Infected"), levels=c("Not infected","Infected"))) %>%
filter(prepost=="Post") %>%
ggplot(aes(x=NMDS1_weighted_unifrac, y=NMDS2_weighted_unifrac)) +
geom_point(aes(bg=time, col=Infection_status), cex=3, pch=21) +
scale_color_manual(values=c("white","red"))+
facet_wrap(~species, nrow=1)
##### Plotting alpha diversity ######
# Fit normal and lognormal distributions to each of the alpha diversity values
# Observed OTUs
x.fit.norm <- seq(min(mf_con$observed_otus)-sd(mf_con$observed_otus)
, max(mf_con$observed_otus)+sd(mf_con$observed_otus)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(mf_con$observed_otus, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(mf_con$observed_otus, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
ggplot(data=mf_con, aes(x=observed_otus)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue")
# Chao1
x.fit.norm <- seq(min(mf_con$chao1)-sd(mf_con$chao1)
, max(mf_con$chao1)+sd(mf_con$chao1)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(mf_con$chao1, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(mf_con$chao1, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
ggplot(data=mf_con, aes(x=chao1)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue")
# Shannon
x.fit.norm <- seq(min(mf_con$shannon)-sd(mf_con$shannon)
, max(mf_con$shannon)+sd(mf_con$shannon)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(mf_con$shannon, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(mf_con$shannon, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
ggplot(data=mf_con, aes(x=shannon)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue")
# Faith's PD
x.fit.norm <- seq(min(mf_con$faith_pd)-sd(mf_con$faith_pd)
, max(mf_con$faith_pd)+sd(mf_con$faith_pd)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(mf_con$faith_pd, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(mf_con$faith_pd, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma distribution
param.gamma <- fitdistr(mf_con$faith_pd, densfun="gamma")
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
ggplot(data=mf_con, aes(x=faith_pd)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green")
# Check to see if turnover is changing with time significantly
options(repr.plot.height=6, repr.plot.width=10)
gg_divtime_con <- mf_con %>%
filter(!is.na(observed_otus)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=observed_otus)) +
geom_line(aes(group=indivID)) +
geom_point(aes(col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_divtime_treat <- mf_treat %>%
filter(!is.na(observed_otus)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=observed_otus)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5)) +
facet_grid(~species)
grid.arrange(gg_divtime_con, gg_divtime_treat, nrow=2)
# Check to see if turnover is changing with time significantly
options(repr.plot.height=6, repr.plot.width=10)
gg_divtime_con <- mf_con %>%
filter(!is.na(chao1)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=chao1)) +
geom_line(aes(group=indivID)) +
geom_point(aes(col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_divtime_treat <- mf_treat %>%
filter(!is.na(chao1)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=chao1)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5)) +
facet_grid(~species)
grid.arrange(gg_divtime_con, gg_divtime_treat, nrow=2)
# Check to see if turnover is changing with time significantly
options(repr.plot.height=6, repr.plot.width=10)
gg_divtime_con <- mf_con %>%
filter(!is.na(shannon)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=shannon)) +
geom_line(aes(group=indivID)) +
geom_point(aes(col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_divtime_treat <- mf_treat %>%
filter(!is.na(shannon)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=shannon)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5)) +
facet_grid(~species)
grid.arrange(gg_divtime_con, gg_divtime_treat, nrow=2)
# Check to see if turnover is changing with time significantly
options(repr.plot.height=6, repr.plot.width=10)
gg_divtime_con <- mf_con %>%
filter(!is.na(faith_pd)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=faith_pd)) +
geom_line(aes(group=indivID)) +
geom_point(aes(col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_divtime_treat <- mf_treat %>%
filter(!is.na(faith_pd)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=faith_pd)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5)) +
facet_grid(~species)
grid.arrange(gg_divtime_con, gg_divtime_treat, nrow=2)
####### Alpha diversity metrics test ######
# Loop through all alpha diversity metrics
for ( a in alpha_metrics ) {
print("-------------")
print(a)
assign(paste0("anova_",a), anova_log_lm_3way(mf=mf_alt_filt_final, dep=paste0(a), indep1="species", indep2="time", indep3="Bd_exposure"))
}
## [1] "-------------"
## [1] "observed_otus"
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(indep1)` instead of `indep1` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(indep2)` instead of `indep2` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(indep3)` instead of `indep3` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## Note: Using an external vector in selections is ambiguous.
## ℹ Use `all_of(dep)` instead of `dep` to silence this message.
## ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This message is displayed once per session.
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.281586333560066, F(4,495) = 1.26794091530175"
## [1] "2-way interaction"
## [1] "species:time p = 1.39010637776659e-06, F(4,499) = 8.42565261946424"
## [1] "species:Bd_exposure p = 0.925624301434885, F(4,499) = 0.222895696278359"
## [1] "time:Bd_exposure p = 0.241741735268746, F(1,499) = 1.37366916620789"
## [1] "1-way main interactions"
## [1] "species p = 5.88947659698197e-18, F(4,495) = 23.6707936679131"
## [1] "time p = 0.336981331871965, F(1,495) = 0.923642088611239"
## [1] "Bd_exposure p = 0.002055338808181, F(1,495) = 9.59912806734385"
## [1] "-------------"
## [1] "chao1"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.448150013097119, F(4,495) = 0.926548227683268"
## [1] "2-way interaction"
## [1] "species:time p = 6.23910894743295e-05, F(4,499) = 6.2741720138224"
## [1] "species:Bd_exposure p = 0.884493187885575, F(4,499) = 0.289981865686314"
## [1] "time:Bd_exposure p = 0.376248405627729, F(1,499) = 0.784324404483501"
## [1] "1-way main interactions"
## [1] "species p = 2.07008549112312e-16, F(4,495) = 21.5255294024166"
## [1] "time p = 0.302035314885298, F(1,495) = 1.06737187902268"
## [1] "Bd_exposure p = 0.00868169653276428, F(1,495) = 6.9413563438679"
## [1] "-------------"
## [1] "shannon"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.356095408056368, F(4,495) = 1.09952504951877"
## [1] "2-way interaction"
## [1] "species:time p = 1.06878066690256e-09, F(4,499) = 12.4987746194204"
## [1] "species:Bd_exposure p = 0.296448730850941, F(4,499) = 1.23158305887371"
## [1] "time:Bd_exposure p = 0.0233406265749638, F(1,499) = 5.17477362046032"
## [1] "1-way main interactions"
## [1] "species p = 1.66208614530443e-06, F(4,495) = 8.3226122585781"
## [1] "time p = 0.000360478097774282, F(1,495) = 12.9019912393834"
## [1] "Bd_exposure p = 0.953121298214154, F(1,495) = 0.00345942895244014"
## [1] "-------------"
## [1] "faith_pd"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.411074583780074, F(4,495) = 0.992646171178883"
## [1] "2-way interaction"
## [1] "species:time p = 6.27209105988508e-08, F(4,499) = 10.1797641084085"
## [1] "species:Bd_exposure p = 0.972632823029778, F(4,499) = 0.127025275816905"
## [1] "time:Bd_exposure p = 0.396669031471287, F(1,499) = 0.719639168926047"
## [1] "1-way main interactions"
## [1] "species p = 4.48226234420104e-21, F(4,495) = 28.0817574637397"
## [1] "time p = 0.192488863650267, F(1,495) = 1.7030209927874"
## [1] "Bd_exposure p = 0.0270163302324145, F(1,495) = 4.9184608014005"
# Now with indivID
for ( a in alpha_metrics ) {
print("-------------")
print(a)
assign(paste0("anova_lmer_",a), anova_log_lmer_3way(mf=mf_alt_filt_final, dep=paste0(a), indep1="species", indep2="time", indep3="Bd_exposure"))
}
## [1] "-------------"
## [1] "observed_otus"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.260364405361499, Chisq(4,4) = 5.2735372283189"
## [1] "2-way interaction"
## [1] "species:time p = 1.34615855229704e-07, Chisq(4,1) = 37.6137389020087"
## [1] "species:Bd_exposure p = 0.988884862092328, Chisq(4,1) = "
## [1] "time:Bd_exposure p = 0.17300373873343, Chisq(1,1) = "
## [1] "1-way main interactions"
## [1] "species p = 1.45749026776731e-15, Chisq(4,4) = 75.6420446919586"
## [1] "time p = 0.299236899405715, Chisq(1,4) = 1.07759185641768"
## [1] "Bd_exposure p = 0.0388996572369675, Chisq(1,4) = 4.26524422254854"
## [1] "-------------"
## [1] "chao1"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.438246296686977, Chisq(4,4) = 3.76846664955745"
## [1] "2-way interaction"
## [1] "species:time p = 2.3200025669242e-05, Chisq(4,1) = 26.667981747491"
## [1] "species:Bd_exposure p = 0.972942629154216, Chisq(4,1) = "
## [1] "time:Bd_exposure p = 0.307658263896612, Chisq(1,1) = "
## [1] "1-way main interactions"
## [1] "species p = 2.41725797326298e-15, Chisq(4,4) = 74.6032703093159"
## [1] "time p = 0.285914974158478, Chisq(1,4) = 1.13875775055916"
## [1] "Bd_exposure p = 0.0448559404661451, Chisq(1,4) = 4.02404756431887"
## [1] "-------------"
## [1] "shannon"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.345940214394255, Chisq(4,4) = 4.47150014462038"
## [1] "2-way interaction"
## [1] "species:time p = 1.86749691571519e-10, Chisq(4,1) = 51.3707346362657"
## [1] "species:Bd_exposure p = 0.630165753494392, Chisq(4,1) = "
## [1] "time:Bd_exposure p = 0.0215048585078992, Chisq(1,1) = "
## [1] "1-way main interactions"
## [1] "species p = 4.96388184796777e-06, Chisq(4,4) = 29.969942983515"
## [1] "time p = 0.000329207280785514, Chisq(1,4) = 12.8964424533697"
## [1] "Bd_exposure p = 0.989371915239734, Chisq(1,4) = 0.000177441656287351"
## [1] "-------------"
## [1] "faith_pd"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.356619432787304, Chisq(4,4) = 4.38322457877476"
## [1] "2-way interaction"
## [1] "species:time p = 4.72248718167167e-09, Chisq(4,1) = 44.640504832635"
## [1] "species:Bd_exposure p = 0.977896246590525, Chisq(4,1) = "
## [1] "time:Bd_exposure p = 0.294061550662723, Chisq(1,1) = "
## [1] "1-way main interactions"
## [1] "species p = 2.51251408986147e-19, Chisq(4,4) = 93.385218661698"
## [1] "time p = 0.198247138653768, Chisq(1,4) = 1.65524072528798"
## [1] "Bd_exposure p = 0.105796766002496, Chisq(1,4) = 2.61591077957358"
# separated controls and Bd_exposures
for ( a in alpha_metrics ) {
print("-------")
print(a)
print("CONTROL")
assign(paste0("anova_", a, "con"), anova_log_lm_2way(mf=mf_con, dep = paste0(a), indep1="species", indep2="time"))
print("Bd_exposure")
assign(paste0("anova_", a, "treat"), anova_log_lm_2way(mf=mf_treat, dep = paste0(a), indep1="species", indep2="time"))
}
## [1] "-------"
## [1] "observed_otus"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.0632314548581625, F(4,206) = 2.26678754684032"
## [1] "1-way main interactions"
## [1] "species p = 7.89810393153774e-23, F(4,210) = 35.8638748376736"
## [1] "time p = 0.0245543160060272, F(1,210) = 5.12869896303142"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 1.25763508404044e-05, F(4,289) = 7.31913638167554"
## [1] "1-way main interactions"
## [1] "species p = 2.11007114225772e-11, F(4,289) = 15.358968133514"
## [1] "time p = 0.416667885543009, F(1,289) = 0.661592533639563"
## [1] "-------"
## [1] "chao1"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.0761962978211546, F(4,206) = 2.14771231845094"
## [1] "1-way main interactions"
## [1] "species p = 1.52854167670683e-19, F(4,210) = 29.6290722527069"
## [1] "time p = 0.0743241802795942, F(1,210) = 3.216813435654"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.000649338599978664, F(4,289) = 5.00449606975323"
## [1] "1-way main interactions"
## [1] "species p = 4.63839186794777e-10, F(4,289) = 13.4419610563413"
## [1] "time p = 0.479089575641407, F(1,289) = 0.502238202181595"
## [1] "-------"
## [1] "shannon"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.00774778327238691, F(4,206) = 3.56707508998551"
## [1] "1-way main interactions"
## [1] "species p = 0.047483544900217, F(4,206) = 2.44802487844419"
## [1] "time p = 0.678091691561947, F(1,206) = 0.172772984469061"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 4.92990449351618e-07, F(4,289) = 9.23291432921763"
## [1] "1-way main interactions"
## [1] "species p = 3.98960258386054e-05, F(4,289) = 6.64107740434005"
## [1] "time p = 0.00266334881616572, F(1,289) = 9.18331528786952"
## [1] "-------"
## [1] "faith_pd"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.0630772965802437, F(4,206) = 2.26834004033497"
## [1] "1-way main interactions"
## [1] "species p = 9.60333834562311e-19, F(4,210) = 28.1798810221786"
## [1] "time p = 0.000438143838525719, F(1,210) = 12.7642818639535"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 5.76584228955603e-07, F(4,289) = 9.13986928014615"
## [1] "1-way main interactions"
## [1] "species p = 1.79225153125794e-14, F(4,289) = 19.8804644476309"
## [1] "time p = 0.250906838194985, F(1,289) = 1.32355070753381"
# separated controls and Bd_exposures, with indivID
for ( a in alpha_metrics ) {
print("-------")
print(a)
print("CONTROL")
assign(paste0("anova_lmer_", a, "con"), anova_log_lmer_2way(mf=mf_con, dep = paste0(a), indep1="species", indep2="time"))
print("Bd_exposure")
assign(paste0("anova_lmer_", a, "treat"), anova_log_lmer_2way(mf=mf_treat, dep = paste0(a), indep1="species", indep2="time"))
}
## [1] "-------"
## [1] "observed_otus"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.0376718481711225, Chisq(4,4) = 10.1692076232745"
## [1] "1-way main interactions"
## [1] "species p = 1.00422742098535e-05, Chisq(4,4) = 28.4642256811462"
## [1] "time p = 0.571501783824005, Chisq(1,4) = 0.320176184402058"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 1.53800096056648e-06, Chisq(4,4) = 32.4635651290777"
## [1] "1-way main interactions"
## [1] "species p = 1.59209139298367e-09, Chisq(4,4) = 46.9101085122364"
## [1] "time p = 0.386624363903054, Chisq(1,4) = 0.749532249045313"
## [1] "-------"
## [1] "chao1"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.0569498347306608, Chisq(4,4) = 9.17162960915955"
## [1] "1-way main interactions"
## [1] "species p = 1.24680748798057e-10, Chisq(4,1) = 52.2099644234552"
## [1] "time p = 0.0558351349394547, Chisq(1,1) = 3.65696859202849"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.000314125074762635, Chisq(4,4) = 21.0176685223584"
## [1] "1-way main interactions"
## [1] "species p = 1.93730319422719e-09, Chisq(4,4) = 46.5008042745184"
## [1] "time p = 0.458202350881037, Chisq(1,4) = 0.550282345475784"
## [1] "-------"
## [1] "shannon"
## [1] "CONTROL"
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## [1] "2-way interaction"
## [1] "species:time p = 0.00648619692882753, Chisq(4,4) = 14.2683003599434"
## [1] "1-way main interactions"
## [1] "species p = 0.044079204415761, Chisq(4,4) = 9.79209951377675"
## [1] "time p = 0.677659036142878, Chisq(1,4) = 0.172772984469304"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 9.24579819615696e-08, Chisq(4,4) = 38.4046168416628"
## [1] "1-way main interactions"
## [1] "species p = 0.000156386373109973, Chisq(4,4) = 22.5407362443649"
## [1] "time p = 0.00250399825900864, Chisq(1,4) = 9.1376696423001"
## [1] "-------"
## [1] "faith_pd"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.0342465348091655, Chisq(4,4) = 10.3969456762524"
## [1] "1-way main interactions"
## [1] "species p = 1.49753400861981e-05, Chisq(4,4) = 27.6080066887845"
## [1] "time p = 0.585225478246589, Chisq(1,4) = 0.297862621223974"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 4.79595390707188e-08, Chisq(4,4) = 39.784576447826"
## [1] "1-way main interactions"
## [1] "species p = 3.95725879493328e-13, Chisq(4,4) = 64.1125315601142"
## [1] "time p = 0.22649520709927, Chisq(1,4) = 1.4627365587335"
#### Plotting dispersion ####
# Bray-curtis
x.fit.norm <- seq(min(mf_con$disper_braycurtis)-sd(mf_con$disper_braycurtis)
, max(mf_con$disper_braycurtis)+sd(mf_con$disper_braycurtis)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(mf_con$disper_braycurtis, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(mf_con$disper_braycurtis, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma
param.gamma <- fitdistr(mf_con$disper_braycurtis, densfun="gamma")
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
ggplot(data=mf_con, aes(x=disper_braycurtis)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green")
# Unweighted Unifrac
x.fit.norm <- seq(min(mf_con$disper_unweighted_unifrac)-sd(mf_con$disper_unweighted_unifrac)
, max(mf_con$disper_unweighted_unifrac)+sd(mf_con$disper_unweighted_unifrac)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(mf_con$disper_unweighted_unifrac, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(mf_con$disper_unweighted_unifrac, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma
param.gamma <- fitdistr(mf_con$disper_unweighted_unifrac, densfun="gamma")
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
ggplot(data=mf_con, aes(x=disper_unweighted_unifrac)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green")
# Weighted Unifrac
x.fit.norm <- seq(min(mf_con$disper_weighted_unifrac)-sd(mf_con$disper_weighted_unifrac)
, max(mf_con$disper_weighted_unifrac)+sd(mf_con$disper_weighted_unifrac)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(mf_con$disper_weighted_unifrac, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(mf_con$disper_weighted_unifrac, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma
param.gamma <- fitdistr(mf_con$disper_weighted_unifrac, densfun="gamma")
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
ggplot(data=mf_con, aes(x=disper_weighted_unifrac)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green")
## Braycurtis
# Check to see if turnover is changing with time significantly
gg_dispertime_con <- mf_con %>%
filter(!is.na(disper_braycurtis)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=disper_braycurtis)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_dispertime_treat <- mf_treat %>%
filter(!is.na(disper_braycurtis)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=disper_braycurtis)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5))+
facet_grid(~species)
grid.arrange(gg_dispertime_con, gg_dispertime_treat, nrow=2)
## Unweighted Unifrac
# Check to see if turnover is changing with time significantly
gg_dispertime_con <- mf_con %>%
filter(!is.na(disper_unweighted_unifrac)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=disper_unweighted_unifrac)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_dispertime_treat <- mf_treat %>%
filter(!is.na(disper_unweighted_unifrac)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=disper_unweighted_unifrac)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5))+
facet_grid(~species)
grid.arrange(gg_dispertime_con, gg_dispertime_treat, nrow=2)
## Weighted Unifrac
# Check to see if turnover is changing with time significantly
gg_dispertime_con <- mf_con %>%
filter(!is.na(disper_weighted_unifrac)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=disper_weighted_unifrac)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_dispertime_treat <- mf_treat %>%
filter(!is.na(disper_weighted_unifrac)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=disper_weighted_unifrac)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5))+
facet_grid(~species)
grid.arrange(gg_dispertime_con, gg_dispertime_treat, nrow=2)
####### Dispersion metrics test ######
# try regular 3-way iterative test
for ( b in beta_metrics ) {
print("--------")
print(b)
assign(paste0("anova_disper_",b), anova_log_lm_3way(mf=mf_alt_filt_final, dep = paste0("disper_",b), indep1="species", indep2="time", indep3="Bd_exposure"))
}
## [1] "--------"
## [1] "braycurtis"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.630085983403839, F(4,495) = 0.645779878589656"
## [1] "2-way interaction"
## [1] "species:time p = 0.00061040058581474, F(4,499) = 4.97645554861041"
## [1] "species:Bd_exposure p = 0.00546418986218119, F(4,499) = 3.7114966033594"
## [1] "time:Bd_exposure p = 0.653052709497064, F(1,499) = 0.202316635538872"
## [1] "1-way main interactions"
## [1] "species p = 0.51624589409124, F(4,495) = 0.814532886716958"
## [1] "time p = 1.01307946358475e-12, F(1,495) = 53.546725992438"
## [1] "Bd_exposure p = 4.17370966378774e-07, F(1,495) = 26.3064904874292"
## [1] "--------"
## [1] "unweighted_unifrac"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.696004962132921, F(4,495) = 0.554253224648829"
## [1] "2-way interaction"
## [1] "species:time p = 0.00201544419101792, F(4,499) = 4.29034421859866"
## [1] "species:Bd_exposure p = 0.0136780837044633, F(4,499) = 3.17170563148351"
## [1] "time:Bd_exposure p = 0.869598017747192, F(1,499) = 0.0269792449123242"
## [1] "1-way main interactions"
## [1] "species p = 0.0261686048491203, F(4,495) = 2.78391660411219"
## [1] "time p = 0.00135166590478791, F(1,495) = 10.3872857670075"
## [1] "Bd_exposure p = 0.456173813676098, F(1,495) = 0.556130897204231"
## [1] "--------"
## [1] "weighted_unifrac"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.0677924740906419, F(4,495) = 2.20126175678126"
## [1] "2-way interaction"
## [1] "species:time p = 0.0765179645743409, F(4,499) = 2.12537538120733"
## [1] "species:Bd_exposure p = 0.0492774317157696, F(4,499) = 2.39876386349435"
## [1] "time:Bd_exposure p = 0.561867134094754, F(1,499) = 0.336932691409228"
## [1] "1-way main interactions"
## [1] "species p = 0.0315775442315638, F(4,495) = 2.67021039068673"
## [1] "time p = 9.44564761991552e-08, F(1,495) = 29.334478294555"
## [1] "Bd_exposure p = 0.00759725611936036, F(1,495) = 7.18381553439865"
# with indivID
for ( b in beta_metrics ) {
print("--------")
print(b)
assign(paste0("anova_disper_lmer_",b), anova_log_lmer_3way(mf=mf_alt_filt_final, dep = paste0("disper_",b), indep1="species", indep2="time", indep3="Bd_exposure"))
}
## [1] "--------"
## [1] "braycurtis"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.631394259403063, Chisq(4,4) = 2.57423352808069"
## [1] "2-way interaction"
## [1] "species:time p = 0.000287491872945836, Chisq(4,1) = 21.2116458272352"
## [1] "species:Bd_exposure p = 0.190378064463001, Chisq(4,1) = "
## [1] "time:Bd_exposure p = 0.713649799931758, Chisq(1,1) = "
## [1] "1-way main interactions"
## [1] "species p = 0.0194699572292092, Chisq(4,4) = 11.7307460647844"
## [1] "time p = 1.90225449669628e-14, Chisq(1,4) = 58.6306952695872"
## [1] "Bd_exposure p = 0.00679488404671445, Chisq(1,4) = 7.32640462412022"
## [1] "--------"
## [1] "unweighted_unifrac"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.687178932425889, Chisq(4,4) = 2.26483265954386"
## [1] "2-way interaction"
## [1] "species:time p = 0.00173632320412176, Chisq(4,1) = 17.2396200369315"
## [1] "species:Bd_exposure p = 0.0230564821895607, Chisq(4,1) = "
## [1] "time:Bd_exposure p = 0.8670093957266, Chisq(1,1) = "
## [1] "1-way main interactions"
## [1] "species p = 0.0287233352713933, Chisq(4,4) = 10.81503457663"
## [1] "time p = 0.00125939107697946, Chisq(1,4) = 10.4011170185069"
## [1] "Bd_exposure p = 0.488714579267407, Chisq(1,4) = 0.479355090818542"
## [1] "--------"
## [1] "weighted_unifrac"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.0512379242080112, Chisq(4,4) = 9.42847180353398"
## [1] "2-way interaction"
## [1] "species:time p = 0.0614523149394471, Chisq(4,1) = 8.98592255296896"
## [1] "species:Bd_exposure p = 0.206464934077034, Chisq(4,1) = "
## [1] "time:Bd_exposure p = 0.539570711046414, Chisq(1,1) = "
## [1] "1-way main interactions"
## [1] "species p = 0.477779332537429, Chisq(4,4) = 3.5006512162275"
## [1] "time p = 2.66482832686956e-08, Chisq(1,4) = 30.9374271094969"
## [1] "Bd_exposure p = 0.1063943575398, Chisq(1,4) = 2.60697768638262"
# Manual for weighted unifrac
Anova(lm(log(disper_weighted_unifrac) ~ species*time*Bd_exposure, data = mf_alt_filt_final), type=3)
## Anova Table (Type III tests)
##
## Response: log(disper_weighted_unifrac)
## Sum Sq Df F value Pr(>F)
## (Intercept) 0.737 1 5.3391 0.021262 *
## species 0.457 4 0.8277 0.507920
## time 0.810 1 5.8633 0.015818 *
## Bd_exposure 1.253 1 9.0718 0.002729 **
## species:time 0.348 4 0.6306 0.640866
## species:Bd_exposure 1.981 4 3.5854 0.006783 **
## time:Bd_exposure 0.553 1 4.0035 0.045953 *
## species:time:Bd_exposure 1.216 4 2.2013 0.067792 .
## Residuals 68.371 495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(lmer(log(disper_weighted_unifrac) ~ species*time*Bd_exposure + (1|indivID), data = mf_alt_filt_final), type=3)
## Analysis of Deviance Table (Type III Wald chisquare tests)
##
## Response: log(disper_weighted_unifrac)
## Chisq Df Pr(>Chisq)
## (Intercept) 4.9504 1 0.02608 *
## species 2.8471 4 0.58373
## time 6.2768 1 0.01223 *
## Bd_exposure 7.9761 1 0.00474 **
## species:time 2.6449 4 0.61890
## species:Bd_exposure 12.5764 4 0.01354 *
## time:Bd_exposure 4.2880 1 0.03838 *
## species:time:Bd_exposure 9.4285 4 0.05124 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
for ( b in beta_metrics ) {
print("-------")
print(b)
print("CONTROL")
assign(paste0("anova_disper_",b,"_con"),anova_log_lm_2way(mf=mf_con, dep = paste0("disper_",b), indep1="species", indep2="time"))
print("Bd_exposure")
assign(paste0("anova_disper_",b,"_treat"),anova_log_lm_2way(mf=mf_treat, dep = paste0("disper_",b), indep1="species", indep2="time"))
}
## [1] "-------"
## [1] "braycurtis"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.0719983796529365, F(4,206) = 2.18398948209259"
## [1] "1-way main interactions"
## [1] "species p = 1.13449473842694e-06, F(4,210) = 8.92031779463094"
## [1] "time p = 7.23788719502915e-08, F(1,210) = 31.1762693787111"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.0101036470931872, F(4,289) = 3.37849289856146"
## [1] "1-way main interactions"
## [1] "species p = 0.882003676298468, F(4,289) = 0.293629658875859"
## [1] "time p = 2.9371611793151e-07, F(1,289) = 27.573017331104"
## [1] "-------"
## [1] "unweighted_unifrac"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.0123981417057312, F(4,206) = 3.28009797529143"
## [1] "1-way main interactions"
## [1] "species p = 0.600602692904124, F(4,206) = 0.688597813546892"
## [1] "time p = 0.0244481935551682, F(1,206) = 5.13785653604999"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.182866618258819, F(4,289) = 1.56802389453821"
## [1] "1-way main interactions"
## [1] "species p = 0.20051548623416, F(4,293) = 1.50563283118577"
## [1] "time p = 0.00060131960463775, F(1,293) = 12.0324273413895"
## [1] "-------"
## [1] "weighted_unifrac"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.00253103395705901, F(4,206) = 4.24536216917606"
## [1] "1-way main interactions"
## [1] "species p = 0.000443482659622865, F(4,206) = 5.29603223917151"
## [1] "time p = 2.13930496506662e-06, F(1,206) = 23.7968206157644"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.677963489077458, F(4,289) = 0.579172122385783"
## [1] "1-way main interactions"
## [1] "species p = 0.0419259618601175, F(4,293) = 2.5118666602411"
## [1] "time p = 0.000243400521952541, F(1,293) = 13.8003837354815"
# with indivID
for ( b in beta_metrics ) {
print("-------")
print(b)
print("CONTROL")
assign(paste0("anova_disper_lmer_",b,"_con"),anova_log_lmer_2way(mf=mf_con, dep = paste0("disper_",b), indep1="species", indep2="time"))
print("Bd_exposure")
assign(paste0("anova_disper_lmer_",b,"_treat"),anova_log_lmer_2way(mf=mf_treat, dep = paste0("disper_",b), indep1="species", indep2="time"))
}
## [1] "-------"
## [1] "braycurtis"
## [1] "CONTROL"
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## [1] "2-way interaction"
## [1] "species:time p = 0.0680489837909238, Chisq(4,4) = 8.7359579283704"
## [1] "1-way main interactions"
## [1] "species p = 3.36516187162997e-07, Chisq(4,1) = 35.6812711785238"
## [1] "time p = 2.35629000590862e-08, Chisq(1,1) = 31.1762693787111"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.00465960674658333, Chisq(4,4) = 15.0201499529926"
## [1] "1-way main interactions"
## [1] "species p = 0.902355864793773, Chisq(4,4) = 1.04849298674853"
## [1] "time p = 1.78553186435506e-08, Chisq(1,4) = 31.7148338763984"
## [1] "-------"
## [1] "unweighted_unifrac"
## [1] "CONTROL"
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## [1] "2-way interaction"
## [1] "species:time p = 0.0107022760911344, Chisq(4,4) = 13.1203919011658"
## [1] "1-way main interactions"
## [1] "species p = 0.599730949338292, Chisq(4,4) = 2.75439125418763"
## [1] "time p = 0.0234095340887611, Chisq(1,4) = 5.13785653605024"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.175009725577779, Chisq(4,4) = 6.34218295799354"
## [1] "1-way main interactions"
## [1] "species p = 0.428278359180316, Chisq(4,1) = 3.83867843706459"
## [1] "time p = 0.000455191997571921, Chisq(1,1) = 12.2908187723035"
## [1] "-------"
## [1] "weighted_unifrac"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.00105800784916662, Chisq(4,4) = 18.3417960704956"
## [1] "1-way main interactions"
## [1] "species p = 0.00104661396408221, Chisq(4,4) = 18.3658110105548"
## [1] "time p = 4.34208218669857e-07, Chisq(1,4) = 25.535986570111"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.659033463602222, Chisq(4,4) = 2.41990313147111"
## [1] "1-way main interactions"
## [1] "species p = 0.179610752446625, Chisq(4,1) = 6.27382629956479"
## [1] "time p = 0.000140140938975688, Chisq(1,1) = 14.5002500553791"
#### Plotting distance ####
# Bray-curtis
# eliminate NAs
beta_values <- mf_con$dist_braycurtis[!is.na(mf_con$dist_braycurtis)]
x.fit.norm <- seq(min(beta_values)-sd(beta_values)
, max(beta_values)+sd(beta_values)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(beta_values, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(beta_values, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma
param.gamma <- fitdistr(beta_values, densfun="gamma")
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
# Fit a beta
param.beta <- fitdistr(beta_values, densfun="beta", start=list(shape1=6,shape2=6))
y.pred.beta <- dbeta(x.fit.norm, shape1=param.beta$estimate[1], shape2 = param.beta$estimate[2])
ggplot(data=mf_con, aes(x=dist_braycurtis)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.beta), aes(x=x, y=y), col="purple")
## Warning: Removed 33 rows containing non-finite values (stat_bin).
# Unweighted Unifrac
# eliminate NAs
beta_values <- mf_con$dist_unweighted_unifrac[!is.na(mf_con$dist_unweighted_unifrac)]
x.fit.norm <- seq(min(beta_values)-sd(beta_values)
, max(beta_values)+sd(beta_values)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(beta_values, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(beta_values, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma
param.gamma <- fitdistr(beta_values, densfun="gamma")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
# Fit a beta
param.beta <- fitdistr(beta_values, densfun="beta", start=list(shape1=6,shape2=6))
y.pred.beta <- dbeta(x.fit.norm, shape1=param.beta$estimate[1], shape2 = param.beta$estimate[2])
ggplot(data=mf_con, aes(x=dist_unweighted_unifrac)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.beta), aes(x=x, y=y), col="purple")
## Warning: Removed 33 rows containing non-finite values (stat_bin).
# Weighted Unifrac
# eliminate NAs
beta_values <- mf_con$dist_weighted_unifrac[!is.na(mf_con$dist_weighted_unifrac)]
x.fit.norm <- seq(min(beta_values)-sd(beta_values)
, max(beta_values)+sd(beta_values)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(beta_values, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(beta_values, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma
param.gamma <- fitdistr(beta_values, densfun="gamma")
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
# Fit a beta
param.beta <- fitdistr(beta_values, densfun="beta", start=list(shape1=6,shape2=6))
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
y.pred.beta <- dbeta(x.fit.norm, shape1=param.beta$estimate[1], shape2 = param.beta$estimate[2])
ggplot(data=mf_con, aes(x=dist_weighted_unifrac)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.beta), aes(x=x, y=y), col="purple")
## Warning: Removed 33 rows containing non-finite values (stat_bin).
## Braycurtis
# Check to see if turnover is changing with time significantly
gg_disttime_con <- mf_con %>%
filter(!is.na(dist_braycurtis)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=dist_braycurtis)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_disttime_treat <- mf_treat %>%
filter(!is.na(dist_braycurtis)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=dist_braycurtis)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5))+
facet_grid(~species)
grid.arrange(gg_disttime_con, gg_disttime_treat, nrow=2)
## Unweighted Unifrac
# Check to see if turnover is changing with time significantly
gg_disttime_con <- mf_con %>%
filter(!is.na(dist_unweighted_unifrac)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=dist_unweighted_unifrac)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_disttime_treat <- mf_treat %>%
filter(!is.na(dist_unweighted_unifrac)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=dist_unweighted_unifrac)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5))+
facet_grid(~species)
grid.arrange(gg_disttime_con, gg_disttime_treat, nrow=2)
## Weighted Unifrac
# Check to see if turnover is changing with time significantly
gg_disttime_con <- mf_con %>%
filter(!is.na(dist_weighted_unifrac)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=dist_weighted_unifrac)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_disttime_treat <- mf_treat %>%
filter(!is.na(dist_weighted_unifrac)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=dist_weighted_unifrac)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5))+
facet_grid(~species)
grid.arrange(gg_disttime_con, gg_disttime_treat, nrow=2)
####### Distance metrics test ######
# No indivID versions for these, because betareg doesn't have an easy mixed effects model
for ( b in beta_metrics) {
print("--------")
print(b)
assign(paste0("anova_dist_",b), anova_betareg_3way(mf=mf_alt_filt_final, dep = paste0("dist_",b), indep1="species", indep2="time", indep3="Bd_exposure"))
}
## [1] "--------"
## [1] "braycurtis"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.858766048125945, Chisq(4) = 1.31534882062691"
## [1] "2-way interaction"
## [1] "species:time p = 0.830179419310584, Chisq(4) = 1.47999316665709"
## [1] "species:Bd_exposure p = 0.354081176527826, Chisq(4) = 4.40401365280203"
## [1] "time:Bd_exposure p = 0.416860211833961, Chisq(1) = 0.659149926026109"
## [1] "1-way main interactions"
## [1] "species p = 1.09910424510229e-09, Chisq(4) = 47.6825539505676"
## [1] "time p = 0.21804219242779, Chisq(1) = 1.5172099980203"
## [1] "Bd_exposure p = 0.849223423270906, Chisq(1) = 0.0361415747511262"
## [1] "--------"
## [1] "unweighted_unifrac"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.988836982189307, Chisq(4) = 0.314827971450356"
## [1] "2-way interaction"
## [1] "species:time p = 0.0796371019133813, Chisq(4) = 8.34780482559164"
## [1] "species:Bd_exposure p = 0.183834777088311, Chisq(4) = 6.21244204283143"
## [1] "time:Bd_exposure p = 0.150918531644647, Chisq(1) = 2.06294216120264"
## [1] "1-way main interactions"
## [1] "species p = 6.47995891357395e-47, Chisq(4) = 222.14385212608"
## [1] "time p = 0.000411449141827737, Chisq(1) = 12.4794689778126"
## [1] "Bd_exposure p = 0.111745812512039, Chisq(1) = 2.52935603235984"
## [1] "--------"
## [1] "weighted_unifrac"
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.979832908024619, Chisq(4) = 0.431324093006434"
## [1] "2-way interaction"
## [1] "species:time p = 0.713844617832534, Chisq(4) = 2.11920680934332"
## [1] "species:Bd_exposure p = 0.325544635998856, Chisq(4) = 4.64635385479841"
## [1] "time:Bd_exposure p = 0.0586354166975787, Chisq(1) = 3.57556732850924"
## [1] "1-way main interactions"
## [1] "species p = 3.48855980285596e-12, Chisq(4) = 59.6187007114131"
## [1] "time p = 0.572802376726634, Chisq(1) = 0.318016043219164"
## [1] "Bd_exposure p = 0.720120535647192, Chisq(1) = 0.128377245175563"
for ( b in beta_metrics ) {
print("-------")
print(b)
print("CONTROL")
assign(paste0("anova_dist_",b,"_con"),anova_betareg_2way(mf=mf_con, dep = paste0("dist_",b), indep1="species", indep2="time"))
print("Bd_exposure")
assign(paste0("anova_dist_",b,"_treat"),anova_betareg_2way(mf=mf_treat, dep = paste0("dist_",b), indep1="species", indep2="time"))
}
## [1] "-------"
## [1] "braycurtis"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.796211181230552, Chisq(4) = 1.66971564653508"
## [1] "1-way main interactions"
## [1] "species p = 1.38764340701417e-06, Chisq(4) = 32.6819517794248"
## [1] "time p = 0.785205319517455, Chisq(1) = 0.0742791291292705"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.863871838162616, Chisq(4) = 1.28525730344115"
## [1] "1-way main interactions"
## [1] "species p = 8.59432493188854e-05, Chisq(4) = 23.841300385566"
## [1] "time p = 0.183649949413768, Chisq(1) = 1.76783735965566"
## [1] "-------"
## [1] "unweighted_unifrac"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.506531077478599, Chisq(4) = 3.31517774558746"
## [1] "1-way main interactions"
## [1] "species p = 4.49048061163884e-28, Chisq(4) = 134.385522574168"
## [1] "time p = 0.000769755801892831, Chisq(1) = 11.3127765634732"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.25073863535415, Chisq(4) = 5.3771757387196"
## [1] "1-way main interactions"
## [1] "species p = 1.54142219094431e-20, Chisq(4) = 99.0835940537224"
## [1] "time p = 0.0840544554964489, Chisq(1) = 2.9847078287114"
## [1] "-------"
## [1] "weighted_unifrac"
## [1] "CONTROL"
## [1] "2-way interaction"
## [1] "species:time p = 0.88872664848568, Chisq(4) = 1.1347103206112"
## [1] "1-way main interactions"
## [1] "species p = 5.01518730839011e-10, Chisq(4) = 49.3164902540102"
## [1] "time p = 0.269640914974441, Chisq(1) = 1.21857283842974"
## [1] "Bd_exposure"
## [1] "2-way interaction"
## [1] "species:time p = 0.839247205952655, Chisq(4) = 1.42838306353759"
## [1] "1-way main interactions"
## [1] "species p = 6.59731123403802e-05, Chisq(4) = 24.4140222438178"
## [1] "time p = 0.161515865431723, Chisq(1) = 1.95997618527918"
#### Plotting inhib prop #####
# Inhibitory proportion
# eliminate NAs
x.fit.norm <- seq(min(mf_con$percInhib)-sd(mf_con$percInhib)
, max(mf_con$percInhib)+sd(mf_con$percInhib)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(mf_con$percInhib, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(mf_con$percInhib, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma
param.gamma <- fitdistr(mf_con$percInhib, densfun="gamma")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
# Fit a beta
param.beta <- fitdistr(mf_con$percInhib, densfun="beta", start=list(shape1=6,shape2=6))
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
y.pred.beta <- dbeta(x.fit.norm, shape1=param.beta$estimate[1], shape2 = param.beta$estimate[2])
ggplot(data=mf_con, aes(x=percInhib)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.beta), aes(x=x, y=y), col="purple")
# Check to see if turnover is changing with time significantly
gg_percInhibtime_con <- mf_con %>%
filter(!is.na(percInhib)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=percInhib)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_percInhibtime_treat <- mf_treat %>%
filter(!is.na(percInhib)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=percInhib)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5))+
facet_grid(~species)
grid.arrange(gg_percInhibtime_con, gg_percInhibtime_treat, nrow=2)
#### Inhib prop metrics test #####
# Check ANOVAs to see if statistical change in percent inhibitory. I think the beta distribution above looks the best.
anova_percInhib <- anova_betareg_3way(mf=mf_alt_filt_final, dep="percInhib", indep1="species", indep2="time", indep3="Bd_exposure")
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.203688453280079, Chisq(4) = 5.93981215805638"
## [1] "2-way interaction"
## [1] "species:time p = 1.11677862256583e-24, Chisq(4) = 118.500196243288"
## [1] "species:Bd_exposure p = 0.0160985724741914, Chisq(4) = 12.1746678239264"
## [1] "time:Bd_exposure p = 0.0715952510467093, Chisq(1) = 3.24606262112476"
## [1] "1-way main interactions"
## [1] "species p = 2.06836863420317e-23, Chisq(4) = 112.561319972827"
## [1] "time p = 0.226725883902658, Chisq(1) = 1.46128431138861"
## [1] "Bd_exposure p = 0.822140226664606, Chisq(1) = 0.0505319421160109"
# Separate control and treat
print("CONTROL")
## [1] "CONTROL"
anova_percInhib_con <- anova_betareg_2way_percInhib(mf=mf_con, dep = "percInhib", indep1="species", indep2="time")
## [1] "2-way interaction"
## [1] "species:time p = 5.63841182540648e-10, Chisq(4) = 49.0727029612388"
## [1] "1-way main interactions"
## [1] "species p = 6.18261617198646e-19, Chisq(4) = 91.5453488168655"
## [1] "time p = 0.0137739292096693, Chisq(1) = 6.06688982957216"
print("Bd_exposure")
## [1] "Bd_exposure"
anova_percInhib_treat <- anova_betareg_2way_percInhib(mf=mf_treat, dep = "percInhib", indep1="species", indep2="time")
## [1] "2-way interaction"
## [1] "species:time p = 3.43433461985765e-15, Chisq(4) = 73.8819702469578"
## [1] "1-way main interactions"
## [1] "species p = 1.09508879526298e-15, Chisq(4) = 76.2288649176893"
## [1] "time p = 0.842929262332258, Chisq(1) = 0.0392626341490796"
##### Plotting inhibitory richness #####
# Inhibitory richness
# eliminate NAs
inhibRich_dat <- mf_con$inhibRich
x.fit.norm <- seq(min(inhibRich_dat)-sd(inhibRich_dat)
, max(inhibRich_dat)+sd(inhibRich_dat)
, length.out = 100)
# Fit normal distribution
param.norm <- fitdistr(inhibRich_dat, densfun="normal")
y.pred.norm <- dnorm(x.fit.norm, mean = param.norm$estimate[1], sd = param.norm$estimate[2])
# Fit a lognormal distribution
param.lnorm <- fitdistr(inhibRich_dat, densfun="lognormal")
y.pred.lnorm <- dlnorm(x.fit.norm, meanlog=param.lnorm$estimate[1], sdlog = param.lnorm$estimate[2])
# Fit a gamma
param.gamma <- fitdistr(inhibRich_dat, densfun="gamma")
## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
y.pred.gamma <- dgamma(x.fit.norm, shape=param.gamma$estimate[1], rate = param.gamma$estimate[2])
# Fit a poisson
param.pois <- fitdistr(inhibRich_dat, densfun="poisson")
y.pred.pois <- dpois(round(x.fit.norm), lambda=param.pois$estimate[1])
ggplot(data=mf_con, aes(x=inhibRich)) +
geom_histogram(aes(y=..density..), bins=20) +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.norm), aes(x=x, y=y), col="red") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.lnorm), aes(x=x, y=y), col="blue") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.gamma), aes(x=x, y=y), col="green") +
geom_line(data=data.frame(x=x.fit.norm, y=y.pred.pois), aes(x=x, y=y), col="purple")
# Check to see if turnover is changing with time significantly
gg_inhibRichtime_con <- mf_con %>%
filter(!is.na(inhibRich)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=inhibRich)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
facet_grid(~species)
gg_inhibRichtime_treat <- mf_treat %>%
filter(!is.na(inhibRich)) %>%
mutate(PABD=factor(PABD)) %>%
ggplot(aes(x=time, y=inhibRich)) +
geom_line(aes(group=indivID)) +
geom_point(aes(group=indivID, col=PABD)) +
scale_color_manual(values=c("blue","red"))+
geom_vline(aes(xintercept=5.5))+
facet_grid(~species)
grid.arrange(gg_inhibRichtime_con, gg_inhibRichtime_treat, nrow=2)
##### Inhib Rich metrics test #####
# Check ANOVAs to see if statistical change in percent inhibitory. I think the beta distribution above looks the best.
anova_inhibRich <- anova_betareg_3way(mf=mf_alt_filt_final, dep="percInhib", indep1="species", indep2="time", indep3="Bd_exposure")
## [1] "3-way interaction"
## [1] "species:time:Bd_exposure p = 0.203688453280079, Chisq(4) = 5.93981215805638"
## [1] "2-way interaction"
## [1] "species:time p = 1.11677862256583e-24, Chisq(4) = 118.500196243288"
## [1] "species:Bd_exposure p = 0.0160985724741914, Chisq(4) = 12.1746678239264"
## [1] "time:Bd_exposure p = 0.0715952510467093, Chisq(1) = 3.24606262112476"
## [1] "1-way main interactions"
## [1] "species p = 2.06836863420317e-23, Chisq(4) = 112.561319972827"
## [1] "time p = 0.226725883902658, Chisq(1) = 1.46128431138861"
## [1] "Bd_exposure p = 0.822140226664606, Chisq(1) = 0.0505319421160109"
print("CONTROL")
## [1] "CONTROL"
anova_inhibRich_con <- anova_pois_2way(mf=mf_con, dep="inhibRich", indep1="species", indep2="time")
## [1] "2-way interaction"
## [1] "species:time p = 1.70092047496161e-05, LR Chisq(4) = 27.3347640458686"
## [1] "1-way main interactions"
## [1] "species p = 3.34024732013679e-05, LR Chisq(4) = 25.8835399629209"
## [1] "time p = 0.00173243232616666, LR Chisq(1) = 9.813418002345"
print("Bd_exposure")
## [1] "Bd_exposure"
anova_inhibRich_treat <- anova_pois_2way(mf=mf_treat, dep="inhibRich", indep1="species", indep2="time")
## [1] "2-way interaction"
## [1] "species:time p = 0.000107499808542882, LR Chisq(4) = 23.355760508068"
## [1] "1-way main interactions"
## [1] "species p = 7.30857417726049e-08, LR Chisq(4) = 38.8991912997543"
## [1] "time p = 0.00316994978254601, LR Chisq(1) = 8.70697391284472"
STOP HERE
# Re-run everything?
RERUN_RICH = FALSE
RERUN_DIST = FALSE
RERUN_DISP = FALSE
RERUN_PERCINHIB = FALSE
RERUN_INHIBRICH = FALSE
##### RICHNESS (observed otus) (I) #######
if ( RERUN_RICH ) {
lmer_log_observed_otus_wtime <- stan_lmer(log(observed_otus) ~ -1 + species*time + (1|indivID), data=mf_con
, prior = normal(0, 10, autoscale = TRUE)
#, family = gaussian(link="log")
, seed = 988735
, adapt_delta = 0.999
)
save(lmer_log_observed_otus_wtime, file="./4_Bayesian_models/lmer_log_observed_otus_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_log_observed_otus_wtime.RData")
}
prior_summary(lmer_log_observed_otus_wtime)
observed_otus_processed_wtime <- process_glmer_withtime(g_lmer = lmer_log_observed_otus_wtime
, dep= "observed_otus"
, name_dep="log_observed_otus"
, transform_raw_func="log"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor=T
, time_inter = T)
observed_otus_processed_wtime$gg_model_Distribution_of_all_values
observed_otus_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
observed_otus_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
observed_otus_processed_wtime$gg_ExpectedDistribution_controls
##### RICHNESS (Chao1 otus) (I) #######
if ( RERUN_RICH ) {
lmer_log_chao1_wtime <- stan_lmer(log(chao1) ~ -1 + species*time + (1|indivID), data=mf_con
, prior = normal(0, 10, autoscale = TRUE)
#, family = gaussian(link="log")
, seed = 5793482
, adapt_delta = 0.999
)
save(lmer_log_chao1_wtime, file="./4_Bayesian_models/lmer_log_chao1_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_log_chao1_wtime.RData")
}
prior_summary(lmer_log_chao1_wtime)
chao1_processed_wtime <- process_glmer_withtime(g_lmer = lmer_log_chao1_wtime
, dep= "chao1"
, name_dep="log_chao1"
, transform_raw_func="log"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T)
chao1_processed_wtime$gg_model_Distribution_of_all_values
chao1_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
chao1_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
chao1_processed_wtime$gg_ExpectedDistribution_controls
##### RICHNESS (Shannon) (I) #######
if ( RERUN_RICH ) {
lmer_shannon_wtime <- stan_lmer(shannon ~ -1 + species*time + (1|indivID), data=mf_con
, prior = normal(0, 10, autoscale = TRUE)
, seed = 5793482
, adapt_delta = 0.999
)
save(lmer_shannon_wtime, file="./4_Bayesian_models/lmer_shannon_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_shannon_wtime.RData")
}
prior_summary(lmer_shannon_wtime)
shannon_processed_wtime <- process_glmer_withtime(g_lmer = lmer_shannon_wtime
, dep= "shannon"
, name_dep="shannon"
, transform_raw_func="None"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_inter=T
, time_factor = T)
shannon_processed_wtime$gg_model_Distribution_of_all_values
shannon_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
shannon_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
shannon_processed_wtime$gg_ExpectedDistribution_controls
##### RICHNESS (Faith's PD) (I) #######
if ( RERUN_RICH ) {
glmer_faith_pd_wtime <- stan_glmer(faith_pd ~ species*time + (1|indivID), data=mf_con
, prior = normal(0, 10, autoscale = TRUE)
, seed = 5793482
, family=Gamma(link="identity")
, adapt_delta = 0.999
)
save(glmer_faith_pd_wtime, file="./4_Bayesian_models/glmer_faith_pd_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_faith_pd_wtime.RData")
}
prior_summary(glmer_faith_pd_wtime)
faith_pd_processed_wtime <- process_glmer_withtime(g_lmer = glmer_faith_pd_wtime
, dep= "faith_pd"
, name_dep="faith_pd"
, transform_raw_func="None"
, intercept_present=T
, fit_distr="Gamma"
, mf_con = mf_con
, mf_treat=mf_treat
, time_inter = T
, time_factor = T
)
faith_pd_processed_wtime$gg_model_Distribution_of_all_values
faith_pd_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
faith_pd_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
faith_pd_processed_wtime$gg_ExpectedDistribution_controls
##### Dispersion (BC) (I) #######
if ( RERUN_DISP ) {
lmer_disper_braycurtis_wtime <- stan_lmer(log(disper_braycurtis) ~ -1 + species*time + (1|indivID)
, data=mf_con
, prior_intercept = normal(location = 0,scale = 5, autoscale = TRUE)
, prior = normal(location=0, scale=5, autoscale=TRUE)
, seed= 29473
)
save(lmer_disper_braycurtis_wtime, file="./4_Bayesian_models/lmer_disper_braycurtis_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_disper_braycurtis_wtime.RData")
}
prior_summary(lmer_disper_braycurtis_wtime)
disper_braycurtis_processed_wtime <- process_glmer_withtime(g_lmer = lmer_disper_braycurtis_wtime
, dep= "disper_braycurtis"
, name_dep="log_disper_braycurtis"
, transform_raw_func="log"
, intercept_present=F
, time_factor = T
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_inter = T
)
disper_braycurtis_processed_wtime$gg_model_Distribution_of_all_values
disper_braycurtis_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
disper_braycurtis_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
disper_braycurtis_processed_wtime$gg_ExpectedDistribution_controls
##### Dispersion (UWU) (I) #######
if ( RERUN_DISP ) {
lmer_disper_unweighted_unifrac_wtime <- stan_lmer(log(disper_unweighted_unifrac) ~ -1 + species*time + (1|indivID)
, data=mf_con
, prior_intercept = normal(location = 0,scale = 5, autoscale = TRUE)
, prior = normal(location=0, scale=5, autoscale=TRUE)
, seed= 29473
)
save(lmer_disper_unweighted_unifrac_wtime, file="./4_Bayesian_models/lmer_disper_unweighted_unifrac_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_disper_unweighted_unifrac_wtime.RData")
}
prior_summary(lmer_disper_unweighted_unifrac_wtime)
disper_unweighted_unifrac_processed_wtime <- process_glmer_withtime(g_lmer = lmer_disper_unweighted_unifrac_wtime
, dep= "disper_unweighted_unifrac"
, name_dep="log_disper_unweighted_unifrac"
, transform_raw_func="log"
, intercept_present=F
, time_factor = T
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_inter = T
)
disper_unweighted_unifrac_processed_wtime$gg_model_Distribution_of_all_values
disper_unweighted_unifrac_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
disper_unweighted_unifrac_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
disper_unweighted_unifrac_processed_wtime$gg_ExpectedDistribution_controls
##### Dispersion (WU) (I) #######
if ( RERUN_DISP ) {
lmer_disper_weighted_unifrac_wtime <- stan_lmer(log(disper_weighted_unifrac) ~ -1 + species*time + (1|indivID)
, data=mf_con
, prior_intercept = normal(location = 0,scale = 5, autoscale = TRUE)
, prior = normal(location=0, scale=5, autoscale=TRUE)
, seed= 29473
)
save(lmer_disper_weighted_unifrac_wtime, file="./4_Bayesian_models/lmer_disper_weighted_unifrac_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_disper_weighted_unifrac_wtime.RData")
}
prior_summary(lmer_disper_weighted_unifrac_wtime)
disper_weighted_unifrac_processed_wtime <- process_glmer_withtime(g_lmer = lmer_disper_weighted_unifrac_wtime
, dep= "disper_weighted_unifrac"
, name_dep="log_disper_weighted_unifrac"
, transform_raw_func="log"
, intercept_present=F
, time_factor = T
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_inter = T
)
disper_weighted_unifrac_processed_wtime$gg_model_Distribution_of_all_values
disper_weighted_unifrac_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
disper_weighted_unifrac_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
disper_weighted_unifrac_processed_wtime$gg_ExpectedDistribution_controls
##### Distance (BC) (I) #######
if ( RERUN_DIST ) {
glmer_dist_braycurtis_wtime <- stan_glmer(dist_braycurtis ~ -1 + species*time + (1|indivID)
, data=mf_con
, family =mgcv::betar
, prior_intercept = normal(location = 0.5,scale = 2.5, autoscale = TRUE)
, prior = normal(location=0.5, scale=2.5, autoscale=TRUE)
, seed= 623445
)
save(glmer_dist_braycurtis_wtime, file="./4_Bayesian_models/glmer_dist_braycurtis_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_dist_braycurtis_wtime.RData")
}
prior_summary(glmer_dist_braycurtis_wtime)
dist_braycurtis_processed_wtime <- process_glmer_withtime(g_lmer = glmer_dist_braycurtis_wtime
, dep= "dist_braycurtis"
, name_dep="dist_braycurtis"
, transform_raw_func="beta"
, intercept_present=F
, fit_distr="Beta"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor=T
, time_inter = T
)
dist_braycurtis_processed_wtime$gg_model_Distribution_of_all_values
## Warning: Removed 33 rows containing missing values (geom_point).
## Warning: Removed 25 rows containing missing values (geom_point).
dist_braycurtis_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
dist_braycurtis_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
dist_braycurtis_processed_wtime$gg_ExpectedDistribution_controls
##### Distance (UWU) (I) #######
if ( RERUN_DIST ) {
glmer_dist_unweighted_unifrac_wtime <- stan_glmer(dist_unweighted_unifrac ~ -1 + species*time + (1|indivID)
, data=mf_con
, family =mgcv::betar
, prior_intercept = normal(location = 0.5,scale = 2.5, autoscale = TRUE)
, prior = normal(location=0.5, scale=2.5, autoscale=TRUE)
, seed= 623445
)
save(glmer_dist_unweighted_unifrac_wtime, file="./4_Bayesian_models/glmer_dist_unweighted_unifrac_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_dist_unweighted_unifrac_wtime.RData")
}
prior_summary(glmer_dist_unweighted_unifrac_wtime)
dist_unweighted_unifrac_processed_wtime <- process_glmer_withtime(g_lmer = glmer_dist_unweighted_unifrac_wtime
, dep= "dist_unweighted_unifrac"
, name_dep="dist_unweighted_unifrac"
, transform_raw_func="beta"
, intercept_present=F
, fit_distr="Beta"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T
)
dist_unweighted_unifrac_processed_wtime$gg_model_Distribution_of_all_values
## Warning: Removed 33 rows containing missing values (geom_point).
## Warning: Removed 25 rows containing missing values (geom_point).
dist_unweighted_unifrac_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
dist_unweighted_unifrac_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
dist_unweighted_unifrac_processed_wtime$gg_ExpectedDistribution_controls
##### Distance (WU) (I) #######
if ( RERUN_DIST ) {
glmer_dist_weighted_unifrac_wtime <- stan_glmer(dist_weighted_unifrac ~ -1 + species*time + (1|indivID)
, data=mf_con
, family =mgcv::betar
, prior_intercept = normal(location = 0.5,scale = 2.5, autoscale = TRUE)
, prior = normal(location=0.5, scale=2.5, autoscale=TRUE)
, seed= 623445
)
save(glmer_dist_weighted_unifrac_wtime, file="./4_Bayesian_models/glmer_dist_weighted_unifrac_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_dist_weighted_unifrac_wtime.RData")
}
prior_summary(glmer_dist_weighted_unifrac_wtime)
dist_weighted_unifrac_processed_wtime <- process_glmer_withtime(g_lmer = glmer_dist_weighted_unifrac_wtime
, dep= "dist_weighted_unifrac"
, name_dep="dist_weighted_unifrac"
, transform_raw_func="beta"
, intercept_present=F
, fit_distr="Beta"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T
)
dist_weighted_unifrac_processed_wtime$gg_model_Distribution_of_all_values
## Warning: Removed 33 rows containing missing values (geom_point).
## Warning: Removed 25 rows containing missing values (geom_point).
dist_weighted_unifrac_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
dist_weighted_unifrac_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
dist_weighted_unifrac_processed_wtime$gg_ExpectedDistribution_controls
##### Percent Inhibitory (I) #######
if ( RERUN_PERCINHIB ) {
glmer_percInhib_wtime <- stan_glmer(percInhib ~ -1 + species*time + (1|indivID)
, data=mf_con
, family =mgcv::betar
, prior_intercept = normal(location = 0.5,scale = 2.5, autoscale = TRUE)
, prior = normal(location=0.5, scale=2.5, autoscale=TRUE)
, seed= 59283
)
save(glmer_percInhib_wtime, file="./4_Bayesian_models/glmer_percInhib_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_percInhib_wtime.RData")
}
prior_summary(glmer_percInhib_wtime)
percInhib_processed_wtime <- process_glmer_withtime(g_lmer = glmer_percInhib_wtime
, dep= "percInhib"
, name_dep="percInhib"
, transform_raw_func="beta"
, intercept_present=F
, fit_distr="Beta"
, mf_con=mf_con
, mf_treat=mf_treat
, time_inter = T
, time_factor=T
)
percInhib_processed_wtime$gg_model_Distribution_of_all_values
percInhib_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
percInhib_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
percInhib_processed_wtime$gg_ExpectedDistribution_controls
##### Inhibitory Richness (I) #######
if ( RERUN_INHIBRICH ) {
glmer_inhibRich_wtime <- stan_glmer(inhibRich ~ species*time + (1|indivID) + (1|SampleID), data=mf_con
, prior = normal(0, 10, autoscale = TRUE)
, family= poisson(link="identity")
, seed = 5423409)
save(glmer_inhibRich_wtime, file="./4_Bayesian_models/glmer_inhibRich_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_inhibRich_wtime.RData")
}
prior_summary(glmer_inhibRich_wtime)
inhibRich_processed_wtime <- process_glmer_withtime(g_lmer = glmer_inhibRich_wtime
, dep= "inhibRich"
, name_dep="inhibRich"
, transform_raw_func="None"
, intercept_present=T
, fit_distr="Poisson"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor=T
, time_inter=T
)
inhibRich_processed_wtime$gg_model_Distribution_of_all_values
inhibRich_processed_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
inhibRich_processed_wtime$gg_ExpectedDistribution_and_Bd_exposed
inhibRich_processed_wtime$gg_ExpectedDistribution_controls
# Combine all p values and save work
all_tabs <- c("observed_otus"
, "chao1"
, "shannon"
, "faith_pd"
, "disper_braycurtis"
, "disper_unweighted_unifrac"
, "disper_weighted_unifrac"
, "dist_braycurtis"
, "dist_unweighted_unifrac"
, "dist_weighted_unifrac"
, "percInhib"
, "inhibRich")
# Change names of columns so they're different
for ( tab in all_tabs) {
temp_tab <- get(paste0(tab,"_processed_wtime"))$all_p
name1 <- paste0("exp_",tab)
name2 <- paste0("p_",tab)
colnames(temp_tab) <- c("indivID",name1, name2,"infect")
assign(paste0("all_p_",tab), temp_tab)
}
# combine all_p's
all_p <- data.frame(indivID=observed_otus_processed_wtime$all_p$indivID, infect=observed_otus_processed_wtime$all_p$infect)
for ( tab in all_tabs ) {
all_p <- get(paste0("all_p_",tab)) %>%
dplyr::select(indivID, paste0("exp_",tab), paste0("p_",tab)) %>%
full_join(all_p, by="indivID")
}
save(all_p, file="./4_Bayesian_models/all_p.RData")
# Combine all p values from con and save work
# Change names of columns so they're different
for ( tab in all_tabs) {
temp_tab <- get(paste0(tab,"_processed_wtime"))$Control_individuals
name1 <- paste0("exp_",tab)
name2 <- paste0("p_",tab)
colnames(temp_tab) <- c("indivID","species","indiv",name1, name2)
assign(paste0("all_p_con",tab), temp_tab)
}
# combine all_p's
all_p_con <- data.frame(indivID=observed_otus_processed_wtime$Control_individuals$indivID)
for ( tab in all_tabs ) {
all_p_con <- get(paste0("all_p_con",tab)) %>%
dplyr::select(indivID, paste0("exp_",tab), paste0("p_",tab)) %>%
full_join(all_p_con, by="indivID")
}
save(all_p_con, file="./4_Bayesian_models/all_p_con.RData")
## Finally, combine all_p and all_p_con
all_p_temp <- all_p %>%
mutate(Bd_exposure="Bd-exposed")
all_p_combined <- all_p_con %>%
mutate(infect=0, Bd_exposure="Control") %>%
rbind(all_p_temp)
save(all_p_combined, file="./4_Bayesian_models/all_p_combined.RData")
# Make a mf with pre-exposure treatment individuals
mf_all_noinfect <- mf_alt_filt_final %>%
filter(!(Bd_exposure=="Bd-exposed"&prepost=="Post"))
##### RICHNESS (observed otus) (II) #######
# There was no effect of time
if ( RERUN_RICH ) {
lmer_log_observed_otus_all_wtime <- stan_lmer(log(observed_otus) ~ -1 + species*time + (1|indivID), data=mf_all_noinfect
, prior = normal(0, 10, autoscale = TRUE)
#, family = gaussian(link="log")
, seed = 298473
, adapt_delta = 0.999
)
save(lmer_log_observed_otus_all_wtime, file="./4_Bayesian_models/lmer_log_observed_otus_all_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_log_observed_otus_all_wtime.RData")
}
prior_summary(lmer_log_observed_otus_all_wtime)
observed_otus_processed_all_wtime <- process_glmer_all_wtime(g_lmer = lmer_log_observed_otus_all_wtime
, dep= "observed_otus"
, name_dep="log_observed_otus"
, transform_raw_func="log"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T)
observed_otus_processed_all_wtime$gg_model_Distribution_of_all_values
observed_otus_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
observed_otus_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
observed_otus_processed_all_wtime$gg_ExpectedDistribution_controls
##### RICHNESS (Chao1) (II) #######
if ( RERUN_RICH ) {
lmer_log_chao1_all_wtime <- stan_lmer(log(chao1) ~ -1 + species*time + (1|indivID), data=mf_all_noinfect
, prior = normal(0, 10, autoscale = TRUE)
#, family = gaussian(link="log")
, seed = 5793482
, adapt_delta = 0.999
)
save(lmer_log_chao1_all_wtime, file="./4_Bayesian_models/lmer_log_chao1_all_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_log_chao1_all_wtime.RData")
}
prior_summary(lmer_log_chao1_all_wtime)
chao1_processed_all_wtime <- process_glmer_all_wtime(g_lmer = lmer_log_chao1_all_wtime
, dep= "chao1"
, name_dep="log_chao1"
, transform_raw_func="log"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T)
chao1_processed_all_wtime$gg_model_Distribution_of_all_values
chao1_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
chao1_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
chao1_processed_all_wtime$gg_ExpectedDistribution_controls
##### DIVERSITY (Shannon) (II) #######
if ( RERUN_RICH ) {
lmer_shannon_all_wtime <- stan_lmer(shannon ~ -1 + species*time + (1|indivID), data=mf_all_noinfect
, prior = normal(0, 10, autoscale = TRUE)
, seed = 5793482
, adapt_delta = 0.999
)
save(lmer_shannon_all_wtime, file="./4_Bayesian_models/lmer_shannon_all_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_shannon_all_wtime.RData")
}
prior_summary(lmer_shannon_all_wtime)
shannon_processed_all_wtime <- process_glmer_all_wtime(g_lmer = lmer_shannon_all_wtime
, dep= "shannon"
, name_dep="shannon"
, transform_raw_func="None"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T)
shannon_processed_all_wtime$gg_model_Distribution_of_all_values
shannon_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
shannon_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
shannon_processed_all_wtime$gg_ExpectedDistribution_controls
##### DIVERSITY (Faith's PD) (II) #######
if ( RERUN_RICH ) {
glmer_faith_pd_all_wtime <- stan_glmer(faith_pd ~ species*time + (1|indivID), data=mf_all_noinfect
, prior = normal(0, 10, autoscale = TRUE)
, seed = 5793482
, family=Gamma(link="identity")
, adapt_delta = 0.999
)
save(glmer_faith_pd_all_wtime, file="./4_Bayesian_models/glmer_faith_pd_all_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_faith_pd_all_wtime.RData")
}
prior_summary(glmer_faith_pd_all_wtime)
faith_pd_processed_all_wtime <- process_glmer_all_wtime(g_lmer = glmer_faith_pd_all_wtime
, dep= "faith_pd"
, name_dep="faith_pd"
, transform_raw_func="None"
, intercept_present=T
, fit_distr="Gamma"
, mf_con = mf_con
, mf_treat=mf_treat
, time_factor=T
, time_inter = T
)
faith_pd_processed_all_wtime$gg_model_Distribution_of_all_values
faith_pd_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
faith_pd_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
faith_pd_processed_all_wtime$gg_ExpectedDistribution_controls
##### DISPERSION (Bray-curtis) (II) #######
if ( RERUN_DISP ) {
lmer_disper_braycurtis_all_wtime <- stan_lmer(log(disper_braycurtis) ~ -1 + species*time + (1|indivID)
, data=mf_all_noinfect
, prior_intercept = normal(location = 0,scale = 5, autoscale = TRUE)
, prior = normal(location=0, scale=5, autoscale=TRUE)
, seed= 29473
)
save(lmer_disper_braycurtis_all_wtime, file="./4_Bayesian_models/lmer_disper_braycurtis_all_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_disper_braycurtis_all_wtime.RData")
}
prior_summary(lmer_disper_braycurtis_all_wtime)
disper_braycurtis_processed_all_wtime <- process_glmer_all_wtime(g_lmer = lmer_disper_braycurtis_all_wtime
, dep= "disper_braycurtis"
, name_dep="log_disper_braycurtis"
, transform_raw_func="log"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_inter = T
, time_factor = T
)
disper_braycurtis_processed_all_wtime$gg_model_Distribution_of_all_values
disper_braycurtis_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
disper_braycurtis_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
disper_braycurtis_processed_all_wtime$gg_ExpectedDistribution_controls
##### DISPERSION (unweighted Unifrac) (II) #######
if ( RERUN_DISP ) {
lmer_disper_unweighted_unifrac_all_wtime <- stan_lmer(log(disper_unweighted_unifrac) ~ -1 + species*time + (1|indivID)
, data=mf_all_noinfect
, prior_intercept = normal(location = 0,scale = 5, autoscale = TRUE)
, prior = normal(location=0, scale=5, autoscale=TRUE)
, seed= 29473
)
save(lmer_disper_unweighted_unifrac_all_wtime, file="./4_Bayesian_models/lmer_disper_unweighted_unifrac_all_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_disper_unweighted_unifrac_all_wtime.RData")
}
prior_summary(lmer_disper_unweighted_unifrac_all_wtime)
disper_unweighted_unifrac_processed_all_wtime <- process_glmer_all_wtime(g_lmer = lmer_disper_unweighted_unifrac_all_wtime
, dep= "disper_unweighted_unifrac"
, name_dep="log_disper_unweighted_unifrac"
, transform_raw_func="log"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T
)
disper_unweighted_unifrac_processed_all_wtime$gg_model_Distribution_of_all_values
disper_unweighted_unifrac_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
disper_unweighted_unifrac_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
disper_unweighted_unifrac_processed_all_wtime$gg_ExpectedDistribution_controls
##### DISPERSION (weighted Unifrac) (II) #######
if ( RERUN_DISP ) {
lmer_disper_weighted_unifrac_all_wtime <- stan_lmer(log(disper_weighted_unifrac) ~ -1 + species*time + (1|indivID)
, data=mf_all_noinfect
, prior_intercept = normal(location = 0,scale = 5, autoscale = TRUE)
, prior = normal(location=0, scale=5, autoscale=TRUE)
, seed= 29473
)
save(lmer_disper_weighted_unifrac_all_wtime, file="./4_Bayesian_models/lmer_disper_weighted_unifrac_all_wtime.RData")
} else {
load("./4_Bayesian_models/lmer_disper_weighted_unifrac_all_wtime.RData")
}
prior_summary(lmer_disper_weighted_unifrac_all_wtime)
disper_weighted_unifrac_processed_all_wtime <- process_glmer_all_wtime(g_lmer = lmer_disper_weighted_unifrac_all_wtime
, dep= "disper_weighted_unifrac"
, name_dep="log_disper_weighted_unifrac"
, transform_raw_func="log"
, intercept_present=F
, fit_distr="normal"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T
)
disper_weighted_unifrac_processed_all_wtime$gg_model_Distribution_of_all_values
disper_weighted_unifrac_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
disper_weighted_unifrac_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
disper_weighted_unifrac_processed_all_wtime$gg_ExpectedDistribution_controls
##### DISTANCE (Bray-curtis) (II) #######
if ( RERUN_DIST ) {
glmer_dist_braycurtis_all_wtime <- stan_glmer(dist_braycurtis ~ -1 + species*time + (1|indivID)
, data=mf_all_noinfect
, family =mgcv::betar
, prior_intercept = normal(location = 0.5,scale = 2.5, autoscale = TRUE)
, prior = normal(location=0.5, scale=2.5, autoscale=TRUE)
, seed= 623445
)
save(glmer_dist_braycurtis_all_wtime, file="./4_Bayesian_models/glmer_dist_braycurtis_all_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_dist_braycurtis_all_wtime.RData")
}
prior_summary(glmer_dist_braycurtis_all_wtime)
dist_braycurtis_processed_all_wtime <- process_glmer_all_wtime(g_lmer = glmer_dist_braycurtis_all_wtime
, dep= "dist_braycurtis"
, name_dep="dist_braycurtis"
, transform_raw_func="beta"
, intercept_present=F
, fit_distr="Beta"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T
)
dist_braycurtis_processed_all_wtime$gg_model_Distribution_of_all_values
## Warning: Removed 58 rows containing missing values (geom_point).
## Warning: Removed 26 rows containing missing values (geom_point).
dist_braycurtis_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 26 rows containing non-finite values (stat_smooth).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 26 rows containing non-finite values (stat_smooth).
## Warning: Removed 26 rows containing missing values (geom_point).
dist_braycurtis_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
## Warning: Removed 26 rows containing missing values (geom_point).
dist_braycurtis_processed_all_wtime$gg_ExpectedDistribution_controls
## Warning: Removed 58 rows containing missing values (geom_point).
##### DISTANCE (unweighted Unifrac) (II) #######
if ( RERUN_DIST ) {
glmer_dist_unweighted_unifrac_all_wtime <- stan_glmer(dist_unweighted_unifrac ~ -1 + species*time + (1|indivID)
, data=mf_all_noinfect
, family =mgcv::betar
, prior_intercept = normal(location = 0.5,scale = 2.5, autoscale = TRUE)
, prior = normal(location=0.5, scale=2.5, autoscale=TRUE)
, seed= 623445
)
save(glmer_dist_unweighted_unifrac_all_wtime, file="./4_Bayesian_models/glmer_dist_unweighted_unifrac_all_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_dist_unweighted_unifrac_all_wtime.RData")
}
prior_summary(glmer_dist_unweighted_unifrac_all_wtime)
dist_unweighted_unifrac_processed_all_wtime <- process_glmer_all_wtime(g_lmer = glmer_dist_unweighted_unifrac_all_wtime
, dep= "dist_unweighted_unifrac"
, name_dep="dist_unweighted_unifrac"
, transform_raw_func="beta"
, intercept_present=F
, fit_distr="Beta"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T
)
dist_unweighted_unifrac_processed_all_wtime$gg_model_Distribution_of_all_values
## Warning: Removed 58 rows containing missing values (geom_point).
## Warning: Removed 26 rows containing missing values (geom_point).
dist_unweighted_unifrac_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 26 rows containing non-finite values (stat_smooth).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 26 rows containing non-finite values (stat_smooth).
## Warning: Removed 26 rows containing missing values (geom_point).
dist_unweighted_unifrac_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
## Warning: Removed 26 rows containing missing values (geom_point).
dist_unweighted_unifrac_processed_all_wtime$gg_ExpectedDistribution_controls
## Warning: Removed 58 rows containing missing values (geom_point).
##### DISTANCE (weighted Unifrac) (II) #######
if ( RERUN_DIST ) {
glmer_dist_weighted_unifrac_all_wtime <- stan_glmer(dist_weighted_unifrac ~ -1 + species*time + (1|indivID)
, data=mf_all_noinfect
, family =mgcv::betar
, prior_intercept = normal(location = 0.5,scale = 2.5, autoscale = TRUE)
, prior = normal(location=0.5, scale=2.5, autoscale=TRUE)
, seed= 623445
)
save(glmer_dist_weighted_unifrac_all_wtime, file="./4_Bayesian_models/glmer_dist_weighted_unifrac_all_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_dist_weighted_unifrac_all_wtime.RData")
}
prior_summary(glmer_dist_weighted_unifrac_all_wtime)
dist_weighted_unifrac_processed_all_wtime <- process_glmer_all_wtime(g_lmer = glmer_dist_weighted_unifrac_all_wtime
, dep= "dist_weighted_unifrac"
, name_dep="dist_weighted_unifrac"
, transform_raw_func="beta"
, intercept_present=F
, fit_distr="Beta"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T
)
dist_weighted_unifrac_processed_all_wtime$gg_model_Distribution_of_all_values
## Warning: Removed 58 rows containing missing values (geom_point).
## Warning: Removed 26 rows containing missing values (geom_point).
dist_weighted_unifrac_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 26 rows containing non-finite values (stat_smooth).
## `geom_smooth()` using formula 'y ~ x'
## Warning: Removed 26 rows containing non-finite values (stat_smooth).
## Warning: Removed 26 rows containing missing values (geom_point).
dist_weighted_unifrac_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
## Warning: Removed 26 rows containing missing values (geom_point).
dist_weighted_unifrac_processed_all_wtime$gg_ExpectedDistribution_controls
## Warning: Removed 58 rows containing missing values (geom_point).
##### PERCENT INHIBITORY (II) #######
if ( RERUN_PERCINHIB ) {
glmer_percInhib_all_wtime <- stan_glmer(percInhib ~ -1 + species*time + (1|indivID)
, data=mf_all_noinfect
, family =mgcv::betar
, prior_intercept = normal(location = 0.5,scale = 2.5, autoscale = TRUE)
, prior = normal(location=0.5, scale=2.5, autoscale=TRUE)
, seed= 59283
)
save(glmer_percInhib_all_wtime, file="./4_Bayesian_models/glmer_percInhib_all_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_percInhib_all_wtime.RData")
}
prior_summary(glmer_percInhib_all_wtime)
percInhib_processed_all_wtime <- process_glmer_all_wtime(g_lmer = glmer_percInhib_all_wtime
, dep= "percInhib"
, name_dep="percInhib"
, transform_raw_func="beta"
, intercept_present=F
, fit_distr="Beta"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor = T
, time_inter = T
)
percInhib_processed_all_wtime$gg_model_Distribution_of_all_values
percInhib_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
percInhib_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
percInhib_processed_all_wtime$gg_ExpectedDistribution_controls
##### INHIBITORY RICHNESS (II) #######
if ( RERUN_INHIBRICH ) {
# Here, I include a random variable of sample because of over-dispersion from poisson
glmer_inhibRich_all_wtime <- stan_glmer(inhibRich ~ species*time + (1|indivID) + (1|SampleID), data=mf_all_noinfect
, prior = normal(0, 10, autoscale = TRUE)
, family= poisson(link="identity")
, seed = 5423409)
save(glmer_inhibRich_all_wtime, file="./4_Bayesian_models/glmer_inhibRich_all_wtime.RData")
} else {
load("./4_Bayesian_models/glmer_inhibRich_all_wtime.RData")
}
prior_summary(glmer_inhibRich_all_wtime)
inhibRich_processed_all_wtime <- process_glmer_all_wtime(g_lmer = glmer_inhibRich_all_wtime
, dep= "inhibRich"
, name_dep="inhibRich"
, transform_raw_func="None"
, intercept_present=T
, fit_distr="Poisson"
, mf_con=mf_con
, mf_treat=mf_treat
, time_factor=T
, time_inter = T
)
# inhibRich_processed_all$gg_model_Distribution_of_all_values
# inhibRich_processed_all$gg_p
# inhibRich_processed_all$gg_ExpectedDistribution_and_Bd_exposed
# inhibRich_processed_all$gg_ExpectedDistribution_controls
inhibRich_processed_all_wtime$gg_model_Distribution_of_all_values
inhibRich_processed_all_wtime$gg_p
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
inhibRich_processed_all_wtime$gg_ExpectedDistribution_and_Bd_exposed
inhibRich_processed_all_wtime$gg_ExpectedDistribution_controls
# Combine all p values and save work
all_tabs <- c("observed_otus"
, "chao1"
, "shannon"
, "faith_pd"
, "disper_braycurtis"
, "disper_unweighted_unifrac"
, "disper_weighted_unifrac"
, "dist_braycurtis"
, "dist_unweighted_unifrac"
, "dist_weighted_unifrac"
, "percInhib"
, "inhibRich")
# Change names of columns so they're different
for ( tab in all_tabs) {
temp_tab <- get(paste0(tab,"_processed_all_wtime"))$all_p
name1 <- paste0("exp_",tab)
name2 <- paste0("p_",tab)
colnames(temp_tab) <- c("indivID","time", name1, name2,"Bd_load")
assign(paste0("all_p_pred_",tab), temp_tab)
}
# combine all_p's
all_p_pred <- data.frame(indivID=observed_otus_processed_all_wtime$all_p$indivID
,time=observed_otus_processed_all_wtime$all_p$time
,infect=observed_otus_processed_all_wtime$all_p$infect)
for ( tab in all_tabs ) {
all_p_pred <- get(paste0("all_p_pred_",tab)) %>%
dplyr::select(indivID,paste0("exp_",tab), paste0("p_",tab), time) %>%
left_join(all_p_pred)
}
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
## Joining, by = c("indivID", "time")
save(all_p_pred, file="./4_Bayesian_models/all_p_pred.RData")
write.table(all_p, file="./4_Bayesian_models/all_p.txt", quote=FALSE, row.names = FALSE
, sep="\t")
write.table(all_p_pred, file="./4_Bayesian_models/all_p_pred.txt", quote=FALSE, row.names = FALSE
, sep="\t")